function [] = runParamRecovery_estimate(runId, ctsTime, replication_idx, Yidx)

%%% Add paths
addpath('../general_funcs');
addpath('../model_funcs');
addpath('../model_funcs/dim_specific_funcs/dim2');
addpath('../model_funcs/dim_specific_funcs/dim3');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');

%%% Get input / output paths
config = getConfig_jointProcess(runId);

indivRunFolders.providerAdoptions  = config.runFolder_providerAdoptions;
indivRunFolders.providerExits      = config.runFolder_providerExits;
indivRunFolders.consumerAdoptions = config.runFolder_consumerAdoptions;

if ctsTime == 2
	myFolder = sprintf('%s/paramRecovery/continuous_time/replication_%03d', config.runFolder, replication_idx);
else if ctsTime == 1
	myFolder = sprintf('%s/paramRecovery/day_level/replication_%03d', config.runFolder, replication_idx);
else
	error('Wrong value of ctsTime.')
end; end;

inputMatFile = sprintf('%s/syntheticData.mat', myFolder);


if Yidx == 1; submodelName = 'providerAdoptions'; end;
if Yidx == 2; submodelName = 'providerExits';  end;
if Yidx == 3; submodelName = 'consumerAdoptions';      end;

outputFolder = createFolderIfNotExist(sprintf('%s/%s', myFolder, submodelName));
outputMatFile = sprintf('%s/results.mat', outputFolder);
outputCSVFile = sprintf('%s/results.csv', outputFolder);


load(inputMatFile, 'data', 'truth');


if Yidx == 1
	%%%%% ESTIMATE PROVIDER-ADOPTIONS MODEL
	data = rmfield(data, 'providerExits');
	data = rmfield(data, 'consumerAdoptions');
	
	disp('========== PROVIDER ADOPTIONS MODEL ==========');
	mydata = data.providerAdoptions;
	mytruth.params = truth.params.providerAdoptions;
	paramNames = get_param_names_providerAdoptions(mydata);
	
	% Initialize params0
	NumParams = mydata.dims.NumParams;
	params0 = mytruth.params;
	
	% Estimate model
	if ctsTime == 2
		[params_hat, covMat, params_ses, modelFits] = estimate_model_by_blocks(3, mydata, 100, params0);
	else
		[params_hat, covMat, params_ses, modelFits] = estimate_model(2, mydata, params0); % works only for discrete-time data
	end
	
	% Display results
	myTable = reportSubmodelResults(paramNames, params_hat, params_ses, modelFits, mytruth);
	disp(myTable);
end

if Yidx == 2
	%%%%% ESTIMATE PROVIDER-EXITS MODEL
	load(inputMatFile, 'data', 'truth');
	data = rmfield(data, 'providerAdoptions');
	data = rmfield(data, 'consumerAdoptions');
	
	disp('========== PROVIDER EXITS MODEL ==========');
	mydata = data.providerExits;
	mytruth.params = truth.params.providerExits;
	paramNames = get_param_names_providerExits(mydata);
	
	% Initialize params0
	NumParams = mydata.dims.NumParams;
	params0 = mytruth.params;
	
	% Estimate model
	if ctsTime == 2
		[params_hat, covMat, params_ses, modelFits] = estimate_model_by_blocks(3, mydata, 100, params0);
	else
		[params_hat, covMat, params_ses, modelFits] = estimate_model(2, mydata, params0); % works only for discrete-time data
	end
	
	% Display results
	myTable = reportSubmodelResults(paramNames, params_hat, params_ses, modelFits, mytruth);
	disp(myTable);
end

if Yidx == 3
	%%%%% ESTIMATE CONSUMER-ADOPTIONS MODEL
	load(inputMatFile, 'data', 'truth');
	data = rmfield(data, 'providerAdoptions');
	data = rmfield(data, 'providerExits');
	
	disp('========== CONSUMER ADOPTIONS MODEL ==========');
	mydata = data.consumerAdoptions;
	mytruth.params = truth.params.consumerAdoptions;
	paramNames = get_param_names_consumerAdoptions(mydata);
	
	% Check: verify that there is at least one transaction per k1, and per k2
	T = size(mydata.Xparts{4}.X,1);
	K1 = size(mydata.Xparts{5}.X,1);
	K2 = size(mydata.Xparts{6}.X,1);
	mydata.Y = reshape(mydata.Y,[T K1*K2]); % T x (K1*K2)
	tmp = full(reshape(sum(mydata.Y,1), [K1 K2])); % K1 x K2
	tmp1 = sum(tmp,2);  % K1 x 1 
	tmp2 = sum(tmp,1)'; % K2 x 1
	FE1 = mydata.Xparts{5}.X_FEs(:,1);
	FE1_NumVals = mydata.Xparts{5}.NumX_FE_vals(1);
	tmp1b = accumarray(FE1, tmp1, [FE1_NumVals 1]); % NumRegions x 1
	disp('Number of consumer adoption per origin region: ');
	disp(tmp1b);
	if ~all(tmp1b>0); error('Some origin region with no transaction!'); end;
	
	FE2 = mydata.Xparts{6}.X_FEs(:,1);
	FE2_NumVals = mydata.Xparts{6}.NumX_FE_vals(1);
	tmp2b = accumarray(FE2, tmp2, [FE2_NumVals 1]); % NumRegions x 1
	disp('Number of consumer adoption per destination region: ');
	disp(tmp2b);
	if ~all(tmp2b>0); error('Some destination region with no transaction!'); end;

	mydata.Y = reshape(mydata.Y,[T*K1*K2 1]); % (T*K1*K2) x 1
	clear tmp tmp1 tmp2 tmp1b tmp2b FE1 FE2 FE1_NumVals FE2_NumVals;
	disp(sprintf('meanY: %d', full(mean(mydata.Y(:)))));
	disp(sprintf('sumY: %d', full(sum(mydata.Y(:)))));
	disp(sprintf('maxY: %d', full(max(mydata.Y(:)))));

	% Initialize params0
	NumParams = mydata.dims.NumParams;
	intercept_idx = find(strcmp(paramNames, 'Beta ttIntercept on NewConsumers'));
	params0 = mytruth.params;
		
	% Estimate model
	if ctsTime == 2
		[params_hat, covMat, params_ses, modelFits] = estimate_model_by_blocks(3, mydata, 100, params0);
	else
		mydata.logLambdaOffset = mydata.logLambdaOffset_dim1 + mydata.logLambdaOffset_dim2; % T x K1 x K2
		mydata.logLambdaOffset = reshape(mydata.logLambdaOffset, [prod(size(mydata.logLambdaOffset)) 1]); % (T*K1*K2) x 1
		[params_hat, covMat, params_ses, modelFits] = estimate_model(2, mydata, params0); % works only for discrete-time data
	end
	% Display results
	myTable = reportSubmodelResults(paramNames, params_hat, params_ses, modelFits, mytruth);
	disp(myTable);
end

%%% Output results to files
save(outputMatFile, 'params_hat', 'params_ses', 'covMat', 'myTable', 'modelFits', '-v7.3');

% Keep 2 decimal points and write to CSV
myTable = formatTableValues(myTable);
writetable(myTable, outputCSVFile,'WriteRowNames',true);

end